# Replication Archive for 
# Alexander Coppock, Emily Ekins and David Kirby (2018), 
# "The Long-lasting Effects of Newspaper Op-Eds on Public Opinion",
# Quarterly Journal of Political Science: Vol. 13: No. 1, pp 59-87.

rm(list = ls())

# Uncomment to install
#install.packages(c("tidyverse", "estimatr"))

library(tidyverse)
library(estimatr)

# Set working directory to replication archive

load("mturk_opeds_cleaned.rdata")
load("elite_opeds_cleaned.rdata")
load("dv_df.rdata")
source("opeds_source.R")


fit_1 <- lm_robust(dv_mj_1_s_w2 ~ Z_distract, mturk_opeds)
fit_2 <- lm_robust(dv_mj_2_s_w2 ~ Z_distract, mturk_opeds)
fit_3 <- lm_robust(dv_mj_3_s_w2 ~ Z_distract, mturk_opeds)
fit_4 <- lm_robust(dv_mj_4_s_w2 ~ Z_distract, mturk_opeds)
fit_5 <- lm_robust(dv_mj_scale_w2 ~ Z_distract, mturk_opeds)
fit_6 <- lm_robust(dv_mj_agree_w2 ~ Z_distract, mturk_opeds)

fit_7 <- lm_robust(dv_mj_1_s_w3 ~ Z_distract, mturk_opeds)
fit_8 <- lm_robust(dv_mj_2_s_w3 ~ Z_distract, mturk_opeds)
fit_9 <- lm_robust(dv_mj_3_s_w3 ~ Z_distract, mturk_opeds)
fit_10 <- lm_robust(dv_mj_4_s_w3 ~ Z_distract, mturk_opeds)
fit_11 <- lm_robust(dv_mj_scale_w3 ~ Z_distract, mturk_opeds)
fit_12 <- lm_robust(dv_mj_agree_w3 ~ Z_distract, mturk_opeds)

fit_13 <- lm_robust(dv_mj_1_s_w2 ~ Z_distract, elite_opeds)
fit_14 <- lm_robust(dv_mj_2_s_w2 ~ Z_distract, elite_opeds)
fit_15 <- lm_robust(dv_mj_3_s_w2 ~ Z_distract, elite_opeds)
fit_16 <- lm_robust(dv_mj_4_s_w2 ~ Z_distract, elite_opeds)
fit_17 <- lm_robust(dv_mj_scale_w2 ~ Z_distract, elite_opeds)
fit_18 <- lm_robust(dv_mj_agree_w2 ~ Z_distract, elite_opeds)

fit_list <- list(fit_1, fit_2, fit_3, fit_4, fit_5, 
                 fit_6, fit_7, fit_8, fit_9, fit_10, 
                 fit_11, fit_12, fit_13, fit_14, fit_15, 
                 fit_16, fit_17, fit_18)

df <-
  bind_rows(
    fit_1 %>% tidy %>% mutate(study = "mturk", wave = "Wave 2", dv = "mj_1"),
    fit_2 %>% tidy %>% mutate(study = "mturk", wave = "Wave 2", dv = "mj_2"),
    fit_3 %>% tidy %>% mutate(study = "mturk", wave = "Wave 2", dv = "mj_3"),
    fit_4 %>% tidy %>% mutate(study = "mturk", wave = "Wave 2", dv = "mj_4"),
    fit_5 %>% tidy %>% mutate(study = "mturk", wave = "Wave 2", dv = "mj_scale"),
    fit_6 %>% tidy %>% mutate(study = "mturk", wave = "Wave 2", dv = "mj_agree"),
    fit_7 %>% tidy %>% mutate(study = "mturk", wave = "Wave 3", dv = "mj_1"),
    fit_8 %>% tidy %>% mutate(study = "mturk", wave = "Wave 3", dv = "mj_2"),
    fit_9 %>% tidy %>% mutate(study = "mturk", wave = "Wave 3", dv = "mj_3"),
    fit_10 %>% tidy %>% mutate(study = "mturk", wave = "Wave 3", dv = "mj_4"),
    fit_11 %>% tidy %>% mutate(study = "mturk", wave = "Wave 3", dv = "mj_scale"),
    fit_12 %>% tidy %>% mutate(study = "mturk", wave = "Wave 3", dv = "mj_agree"),
    fit_13 %>% tidy %>% mutate(study = "elite", wave = "Wave 2", dv = "mj_1"),
    fit_14 %>% tidy %>% mutate(study = "elite", wave = "Wave 2", dv = "mj_2"),
    fit_15 %>% tidy %>% mutate(study = "elite", wave = "Wave 2", dv = "mj_3"),
    fit_16 %>% tidy %>% mutate(study = "elite", wave = "Wave 2", dv = "mj_4"),
    fit_17 %>% tidy %>% mutate(study = "elite", wave = "Wave 2", dv = "mj_scale"),
    fit_18 %>% tidy %>% mutate(study = "elite", wave = "Wave 2", dv = "mj_agree")
  )

gg_df <-
  df %>%
  filter(term != "(Intercept)") %>%
  mutate(
    study = factor(
      study,
      levels = c("mturk", "elite"),
      labels = c("Mechanical Turk", "Elites")
    ),
    label = factor(
      dv,
      levels = c("mj_1", "mj_2", "mj_3", "mj_4", "mj_scale", "mj_agree"),
      labels = c(
        "Favor legalization",
        "Should not prosecute",
        "Less dangerous than heroin",
        "Predict decreased use if legalized",
        "Composite scale",
        "Agreement (Binary)"
      )
    ),
    label = factor(label, levels = rev(levels(label)))
  )


g <- 
ggplot(gg_df, aes(x = estimate, y = label)) +
  geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0) +
  geom_vline(xintercept = 0, linetype = "dotted") +
  theme_bw() +
  facet_grid(study ~ wave) +
  xlab("Standardized Treatment Effect Estimates") +
  theme(strip.background = element_blank(),
        axis.title.y = element_blank())

g
